######## INFO ########

# PROJECT
# Replication files for: Economic Impact of the South China Sea Dispute in China-Philippines
#                        Relations from 2012 to 2016: A DM-LFM analysis  

# Authors: Dianyi Yang
# R Script
# Purpose: This script produces the results for the MCMC diagnostics.
# Created: 27 Oct 2023
# Updated: 29 Nov 2023
# Inputs: custom_functions/DM-LFM.r
#         stored/dmlfm_import.rdata
#         stored/dmlfm.rdata
# Outputs: Table 8: Summary table for MCMC.
#          Table 9: Quantiles of MCMC simulations.
#          output/MCMCtraceplot.png
#          output/MCMCtraceplot.pdf
#          output/MCMCtraceplot_import.png
#          output/MCMCtraceplot_import.pdf
#          output/MCMCdensplot.png
#          output/MCMCdensplot.pdf
#          output/MCMCdensplot_import.png
#          output/MCMCdensplot_import.pdf
#          output/MCMCgeweke.png
#          output/MCMCgeweke.pdf
#          output/MCMCgeweke_import.png
#          output/MCMCgeweke_import.pdf

######## SETUP ########
rm(list = ls()) # clear workspace

need <- c('coda','bpCausal','tidyverse','kableExtra') # list packages needed
have <- need %in% rownames(installed.packages()) # checks packages you have
if(any(!have)) install.packages(need[!have]) # install missing packages
#devtools::install_github("synth-inference/synthdid") #manually install sdid if the codes above fail to do so automatically
invisible(lapply(need, library, character.only=T)) # load needed packages

setwd(gsub('/code','',dirname(rstudioapi::getSourceEditorContext()$path))) #set wd to root directory
list.files()
source('custom_functions/DM-LFM.r')

######## SUMMARY STATS ########
load('stored/dmlfm_import.rdata')
result.use_import <- effSummary(out) #extract effect information
load('stored/dmlfm.rdata')
result.use <- effSummary(out) #extract effect information

mc_att <- mcmc(result.use$eff_avg_i)
mc_att_import <- mcmc(result.use_import$eff_avg_i)
sum <- rbind(t(summary(mc_att)$statistics), t(summary(mc_att_import)$statistics))
sum <- cbind('Trade'=c('Exports','Imports'),as.data.frame(sum))
kbl(sum, position='h', centering=T, digits=6,
    caption = 'Summary table for MCMC.') #Table 8
######## QUANTILE STATS #######
quant <- rbind(t(summary(mc_att)$quantiles), t(summary(mc_att_import)$quantiles))
quant <- cbind('Trade'=c('Exports','Imports'),as.data.frame(quant))
kbl(quant, position='h', centering=T, digits=6,
    caption = 'Quantiles.')#Table 9
######## TRACE PLOTS #########
png('output/MCMCtraceplot.png', units='in', width=8, height=4, res=300)
traceplot(mc_att) #Figure 18
dev.off()

pdf('output/MCMCtraceplot.pdf', width=8, height=4)
traceplot(mc_att) #Figure 18
dev.off()

png('output/MCMCtraceplot_import.png', units='in', width=8, height=4, res=300)
traceplot(mc_att_import) #Figure 19
dev.off()

pdf('output/MCMCtraceplot_import.pdf', width=8, height=4)
traceplot(mc_att_import) #Figure 19
dev.off()
######## DENSITY PLOTS #########
png('output/MCMCdensplot.png', units='in', width=8, height=4, res=300)
densplot(mc_att) #Figure 20
dev.off()

pdf('output/MCMCdensplot.pdf', width=8, height=4)
densplot(mc_att) #Figure 20
dev.off()

png('output/MCMCdensplot_import.png', units='in', width=8, height=4, res=300)
densplot(mc_att_import) #Figure 21
dev.off()

pdf('output/MCMCdensplot_import.pdf', width=8, height=4)
densplot(mc_att_import) #Figure 21
dev.off()
######## GEWEKE PLOTS #########
geweke.diag(mc_att)
pnorm(geweke.diag(mc_att)$z, lower.tail=T)*2

geweke.diag(mc_att_import)
pnorm(geweke.diag(mc_att_import)$z, lower.tail=F)*2

png('output/MCMCgeweke.png', units='in', width=8, height=4, res=300)
geweke.plot(mc_att, nbins=21) #Figure 22
dev.off()

pdf('output/MCMCgeweke.pdf', width=8, height=4)
geweke.plot(mc_att, nbins=21) #Figure 22
dev.off()

png('output/MCMCgeweke_import.png', units='in', width=8, height=4, res=300)
geweke.plot(mc_att_import, nbins=21) #Figure 23
dev.off()

pdf('output/MCMCgeweke_import.pdf', width=8, height=4)
geweke.plot(mc_att_import, nbins=21) #Figure 23
dev.off()